Temperature scaling in a dense vibro-fluidised granular material 
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The leading order "temperature" of a dense two dimensional gran- 
ular material fluidised by external vibrations is determined. The grain 
interactions are characterised by inelastic collisions, but the coeffi- 
cient of restitution is considered to be close to 1, so that the dissi- 
pation of energy during a collision is small compared to the average 
energy of a particle. An asymptotic solution is obtained where the 
particles are considered to be elastic in the leading approximation. 
The velocity distribution is a Maxwell-Boltzmann distribution in the 
leading approximation. The density profile is determined by solving 
the momentum balance equation in the vertical direction, where the 
relation between the pressure and density is provided by the virial 
equation of state. The temperature is determined by relating the 
source of energy due to the vibrating surface and the energy dis- 
sipation due to inelastic collisions. The predictions of the present 
analysis show good agreement with simulation results at higher den- 
sities where theories for a dilute vibrated granular material, with the 
pressure-density relation provided by the ideal gas law, are in error. 



I. INTRODUCTION 

Recent developments in the physics of granular matter ^ 
have illustrated that the dissipative nature of the interactions 
between grains can result in a variety of different phenomena. 
Of particular interest in recent years has been the dynamics 
of vibrated granular materials which exhibit stationary 
states as well as waves and complex patterns. In order to de- 
scribe these diverse states of the material, it is necessary to 
derive macroscopic descriptions by averaging over the micro- 
scopic details of the motion and interactions between individ- 
ual grains. This goal has proved elusive, however, because a 
vibrated granular material is a driven dissipative system, and 
the interactions between the particles are characterised by a 
loss of energy due to inelastic collisions. The statistical me- 
chanics framework developed for equilibrium or near equilib- 
rium systems cannot be used in this case. Consequently, phe- 
nomenological models [Q-^ have been used to describe the 
dynamics of granular materials. The kinetic theories devel- 
oped for granular flows usually assume that the system 
is close to "equilibrium" and the velocity distribution function 
is close to the Maxwell-Boltzmann distribution. 

Experimental studies and computer simulations have re- 
ported the presence of a uniformly fluidised state in a vibrated 
bed of granular material. Luding, Herrmann and Blumen 
carried out 'Event Driven' (ED) simulations of a two dimen- 
sional system of inelastic disks in a gravitational field vibrated 
from below, and obtained scaling laws for the density varia- 
tions in the bed. An experimental study of a vibrated fluidised 
bed was carried out by Warr, Huntley and Jacques . Their 
experimental set up consisted of steel spheres confined be- 



tween two glass plates that are separated by a distance slightly 
larger than the diameter of the spheres. The particles were flu- 
idised by a vibrating surface at the bottom of the bed, and the 
statistics of the velocity distribution of the particles were ob- 
tained using visualisation techniques. Profiles for the density 
and the mean square velocity were obtained, and the parti- 
cle velocity distributions were also determined at certain po- 
sitions in the bed. Both of these studies reported that there is 
an exponential dependence of the density on the height near 
the top of the bed, similar to the Boltzmann distribution for the 
density of a gas in a gravitational field. However, the depen- 
dence of the density deviates from the exponential behaviour 
near the bottom. The dependence of the mean square veloc- 
ity on the vibration frequency and amplitude were found to be 
different in the two studies. 

A theoretical calculation of the distribution function in a 
vibro-fluidised bed was carried out by Kumaran The 
limit of low dissipation, where the coefficient of restitution e 
is close to 1 was considered. In this limit, the mean square 
velocity of the particles is large compared to the mean square 
of the velocity of the vibrating surface, and the dissipation 
of energy during a binary collision is small compared to the 
energy of a particle. A perturbation approximation is used, 
where the energy dissipation is neglected in the leading order 
approximation, and the system resembles a gas at equilibrium 
in a gravitational field. The velocity distribution function is 
a Maxwell-Boltzmann distribution, and the density decreases 
exponentially from the vibrating surface. The first order cor- 
rection to the distribution due to dissipative effects was cal- 
culated using the moment expansion method, and the results 
were found to be in qualitative agreement with the experi- 
ments of Warr et. al. 

The theoretical predictions [|^,p^ were compared with pre- 
vious experimental and simulation studies by McNamara and 
Luding ^V^]- They found that the theory was in good agree- 
ment with experiments for dilute beds, where the area fraction 
of the particles is low, but there were systematic deviations 
from the theoretical predictions as the area fraction increases. 
This is to be expected, since the analysis assumed that the den- 
sity is small and the pair distribution function was set equal to 
1 and therefore the pressure is related to the density by the 
ideal gas law. These assumptions become inaccurate as the 
area fraction of the bed increases. An approximate method 
for including the correction to the pair distribution function 
was suggested by Huntley [p^. 

In the present analysis, the correction to the low density the- 
ory of Kumaran is determined for a vibro-fluidised bed 
where the coefficient of restitution is close to 1. An asymp- 
totic analysis is used, where the dissipation is neglected in the 
leading approximation. The leading order density and veloc- 
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ity profiles are determined using the momentum balance equa- 
tion in the vertical direction. In contrast to the earlier theory 
lP,p^, the virial equation of state for a non-ideal two dimen- 
sional gas is used to determine the leading order density pro- 
file. The density profile differs from the Boltzmann distribu- 
tion, but the velocity distribution function is still a Maxwell- 
Boltzmann distribution. The leading order temperature is de- 
termined by a balance between the source and dissipation of 
energy as before. The complete equilibrium pair distribution 
function is used to determine the rate of dissipation of energy 
due to inelastic collisions. The results are compared with hard 
sphere MD simulations, and also with earlier theoretical and 
simulation studies. 



and V is the area fraction corresponding to p. If the coefficient 
of restitution is set equal to 1 in the leading approximation, the 
equation for the pressure reduces to the standard virial equa- 
tion of state 



(5) 



The resulting equation from Eq. for the density profile is a 
first order ordinary differential equation, which can be solved 
using the mass conservation condition 



(6) 



11. ANALYSIS 

The system consists of a bed of circular disks (of diameter 
cr) in a gravitational field driven by a vibrating surface. The 
vibrating surface has a periodic amplitude function but no as- 
sumption is made regarding the form of the function. There is 
a source of energy at the vibrating surface due to particle col- 
lisions with the surface, and the dissipation is due to inelastic 
collisions. A balance between the two determines the "tem- 
perature", which is the mean square velocity of the particles. 

The limit of low dissipation, where the coefficient of resti- 
tution e is close to 1, is considered. In this limit, it can be 
shown that the mean square velocity of the particles is large 
compared to the mean square velocity of the vibrating sur- 
face. An asymptotic expansion in the parameter e = [/q /Tq is 
used [^. If the source and dissipation of energy are neglected 
in the leading approximation, the system resembles a gas of 
hard disks at equilibrium in a gravitational field. The velocity 
distribution function is a Maxwell-Boltzmann distribution at 
equilibrium 
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where Tq is the leading order temperature. The density profile 
is determined by solving the momentum balance equation in 
the vertical direction 



- 99 = 0, 
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where p is the pressure, p is the density (number of particles 
per area) and g is the acceleration due to gravity. For a gas at 
equilibrium, the pressure is related to the density by the virial 
equation of state, which in the case of inelastic circular disks 
is 



1 



+ il + e)go{v)v 



(3) 



where .go(i') is the pair distribution function at contact, which 
for circular disks is given by [ ]l3[ ] 
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where N is the number of particles per unit width of the bed. 
Note that the leading order temperature Tq is still unknown 
at this stage. This is determined using a balance between the 
source and dissipation of energy. The source of energy due 
to particle collisions with the vibrating surface is determined 
using an equilibrium average over the increase in energy due 
to particle collisions with the vibrating surface (HIK 
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Here {U'^) represents the mean square velocity of the vibrat- 
ing surface. The rate of dissipation of energy per unit width is 
calculated by averaging over the energy loss over all the colli- 
sions between particles and integrating over the height of the 
bed m 
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Note that the go appearing in Sq and Dq is the Enskog factor 
which accounts for the increase in the frequency of collision 
for hard disks at high densities. The temperature Tq can now 
be determined from the relation 



Sq = Do 
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An analytical solution to the density variation Eq. (||) can be 
determined in the low density limit using the equation of state 
for an ideal gas for the pressure [||]. 
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where the leading order temperature is given by, 
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In the low-density limit the density decays exponentially from 
the bottom of the bed. At higher densities the solution to the 
density variation is no longer exponential throughout, and has 
to be obtained numerically by an iterative scheme. However, 
at large distances from the bottom, the bed is dilute and the 
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ideal gas law holds good, hence the decay is exponential, even 
though near the bottom it is not. This gives a convenient start- 
ing point for the numerical integration from a finite height, 
above which we assume the asymptotic solution (z oo) to 
be given by an exponential decay known to within two unde- 
termined constants. A value for the density and the tempera- 
ture is assumed at this height and the integration is carried out 
up to the vibrating plate (z ~ 0). The complete density profile 
is obtained by combining the numerical and the asymptotic so- 
lutions. If the conditions Eq. (||) and Eq. (H) are not satisfied 
after one such integration, a new value is determined for the 
density and temperature using the Newton-Raphson method, 
and the iteration is repeated till convergence. In cases where 
the convergence is poor, the solution is obtained by continuing 
a low density solution in a parameter such as Na or Uq. 

Viscous dissipation: The above analysis can be easily ex- 
tended to the case of dissipation purely due to viscous drag. 
The expression for the source of energy remains the same as 
given by Eq. (0). A drag law given by — —fiUi is assumed. 
The total leading order rate of dissipation per unit width will 
then be 

oo 

Ddo — M Jdz p Jdu F{u) u • u 
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Unlike Eq. (g|), the leading order dissipation is the same for 
the low density and the high density cases. Nevertheless, the 
density profile has to be obtained numerically in the manner 
outhned above, with Eq. (|l2|) substituted for Eq. (^j) in Eq. (^. 



in. SIMULATION AND RESULTS 

The hard sphere molecular dynamics (MD), also known as 
event driven (ED) method [|[ is used for the simulations of the 
vibro-fluidised bed. Periodic boundary conditions are used in 
the horizontal direction and the vibrating surface at the bottom 
has a sawtooth form for the amplitude function. The simula- 
tions are carried out only for the case of inelastic collisions, 
since the viscous drag requires a different treatment than the 
ED method. 

The density profiles obtained using the present analysis, as 
well as the earlier low density approximations of Kumaran 
[^], are compared with the simulation results in Figs. |l| and ^ 
It is seen that the density profiles of the present analysis are 
in good agreement with the simulation results even when the 
density near the bottom of the bed becomes large, while the 
profiles from the low density approximation have significant 
errors. Fig. ^ shows the nature of the density profile in the 
high density limit in the case of dissipation due to viscous 
drag. Here too the present analysis gives reasonable values 
for packing fraction near the bottom, while the low density 
theory predicts physically incorrect values. 

In a recent work, McNamara and Luding [ pi] ] reported the 
scaling of dissipation with the center of mass obtained from 
simulations. The results agreed with the low density theory of 
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FIG. 1. Exponential decay of packing fraction (u) with a nor- 
malised height (z/a) at low densities. The predictions of the present 
analysis (solid line) and the low density theory (dotted line) of ||8|| is 
compared with simulation (points). Both the predictions are nearly 
identical. Here, e — 0.3, Na — 3, g — 1, and Uo = 6. 



JlO| ] but a systematic deviation was observed at high densities 
in all the cases. This deviation is captured in the present anal- 
ysis. The leading order dissipation at low densities in the bed 
is given by [M] 
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In [ pi] ] the total dissipation obtained from the simulation was 
normalised by a factor taken out from this leading order dissi- 
pation and a non dimensional number was defined as 
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The scaling of this factor with the height of the center of mass 
(h) above the position at rest (ho) was studied. This factor 
was found out for different parameter sets by varying the bot- 
tom wall velocity Uq over several decades such that the bed 
is taken from a densely packed regime to a very low density 
regime. They chose a central data set and varied the parame- 
ters one at a time. It was found that in all the cases considered, 
the scaling relation collapsed to a single curve. The central pa- 
rameter set has the following values = 3.2, a — 1, g — 1, 
e = 0.95. 

The present analysis is valid when e = Uq/Tq <C 1 
and when the frequency of particle-particle collision is much 
greater than the frequency of particle-wall collisions. It can 
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FIG. 2. Deviation of the density profile from the exponential de- 
cay at high densities in the case of dissipation due to inelastic colli- 
sions. The simulation result (points) is captured by the present anal- 
ysis (solid line) which is lower than the exponential decay (dotted 
line) of the low density theory of near the bottom of the bed. 
Here e = 0.3, TVcr = 3, g = 1, and f/o = 1. 
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FIG. 3. Deviation of the density profile from the exponential de- 
cay at high densities in the case of dissipation due to viscous drag. 
The present analysis (solid line) gives physically plausible values for 
the packing fraction near the bottom, while the low density theory 
(dotted line) of |^] predicts values higher than the maximum closed 
packing. Here e = 0.2, Na = 20, g = 20, = 0.1, and Uo = 5. 



be shown that in the leading order the ratio of the frequency 
of particle-particle collisions to the frequency particle-wall 
collisions is \/2'k Na. Hence the present analysis will hold 
good when Na ^1/ V2n. The central set corresponds to 
e = 0.35, Na = 3.2 and therefore we expect the present anal- 
ysis to hold good for this case. Most of the parameter sets 
used in Jill ] also fall within the limits of the theory derived 
here. 

Fig. ^ shows the theoretical predictions of the total dissipa- 
tion for the different cases reported in Fig. 2 in [Mn^. It is com- 
pared with the results of two simulations in FigT^ It is seen 
that the present analysis correctly predicts the lowering of the 
coefficient Cpp at high densities. This reduction in the dissi- 
pation from the constant value at low densities is the net result 
of two opposing factors: (i) decrease in the density from the 
exponential behaviour near the vibrating bottom (see Fig. 
hence reducing the total value of the dissipation, and (ii) in- 
crease in frequency of collisions at high densities, increasing 
the dissipation. 

It is also seen that not all the theoretical predictions collapse 
on to a curve as is the case with the data from the simulation. 
In two of the cases the theory does not agree with the simula- 
tions because (i) in one the value of the perturbation parameter 
is high (e — 1.73) and the leading order theory is valid only 
for low e, and (ii) in the other case the value of Na — 0.65 is 
low. 



In Fig. ^ the apparent mismatch with 'e-' is not a discrep- 
ancy with the model, but has got to do with the formula cho- 
sen used in [ [Tl] ] for the normalisation of the dissipation factor 
Cpp. They had chosen to normalise the dissipation by a fac- 
tor (1 — e). While this might have given a better fit for high 
densities (low center of mass), the correct factor for very low 
densities is (1 — e^) as given by Eq. (pj|). The difference is 
more pronounced in the case of e ^ 1, which, here, has a 
value e = 0.75. A close inspection of the curves 'e-' in Fig. ^ 
and Fig. ^ show that the theory and simulation do indeed agree 
with each other. 

We also note here that the data taken from the reported sim- 
ulation [ [Til ] is for asymmetric sawtooth vibration, whereas our 
simulation is for the symmetric sawtooth. Both these give 
similar results for the scaling of Cpp. Also the theoretical pre- 
dictions for the symmetric and the asymmetric sawtooth are 
identical, indicating that the form of the bottom wall vibration 
does not affect the scaling of the dissipation with the center of 
mass. 



IV. CONCLUSION 

In summary, a theory to describe the state of a vibro- 
fluidised bed in the dense limit was derived. This is different 
from the earlier theory of Kumaran [^,|l^, which is valid for 
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FIG. 4. Theoretical scaling of the normalised dissipation (Cpp) 
against the center of mass (h) above the position at rest {ho) for the 
different cases reported in Jl l[]. All except two — (N+) with e = 1.73 
and (N-) with No — 0.65 collapse on to a single curve in the lin- 
ear region. The parameters indicated correspond to A'^ = 16 (N+), 
N = 0.65 (N-), g = 25 (g+), g = 0.04 (g-), e = 0.99 (e+), 
e = 0.75 (e-), rest of the parameters being same as the one in the 
central set, which has the following values N — 3.2, a = 1, g = 1, 
e = 0.95. 
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FIG. 5. Scaling of the normalised dissipation with the center of 
mass: Predictions from the present analysis is compared with the 
results from our simulations and the reported results in [|ll||. The 
linear portion of all the curves from theory, except two, fall on the 
solid line denoted as 'Theory'. The two exceptions are also shown. A 
set of points correspond to the simulation data with parameter values 
iV = 16 (N+), TV = 0.65 (N-), g = 25 (g+), g = 0.04 (g-), 
e — 0.99 (e+), e = 0.75 (e-); rest of the parameters in a set being 
the same as the one in the central set, which has the following values 
N = 3.2, a = l,g = l,e = 0.95. 



low densities where the ideal gas equation was used and the 
pair distribution function was set equal to 1. We have made 
use of the virial equation of state to obtain a correction to the 
exponential density profile obtained in low densities and the 
pair distribution function is used to calculate the increased fre- 
quency of collisions in the source and the dissipation of en- 
ergy. The theoretical predictions of density and temperature 
were compared with the results obtained from MD simula- 
tion of two dimensional disks. The theory correctly predicts 
the lowering of the density from the exponential value at high 
densities near the bottom. The theory also predicts the scaling 
relations of the total dissipation in the bed reported in Jill]. 
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